source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')
library(GEOquery)

# BGI -----------
sobj_bgi <- read_rds('DE_cells/data/BGI_cleaned.rds')

sobj_bgi$Stage |> unique()

sobj_bgi <- sobj_bgi |>
  mutate(group = str_extract(batch, 'Ctrl|COV|Flu'))

sobj_bgi |> VlnPlot('TRPM2', group.by = 'group', pt.size = 0)

sobj_bgi <- sobj_bgi |>
  UpdateSeuratObject()

m2_fluvhc <- sobj_bgi |>
  FindMarkers(group.by = 'group', ident.1 = 'Flu', ident.2 = 'Ctrl',
              min.pct = 0.001)

m2_fluvhc |>
  as_tibble(rownames = 'gene') |>
  mutate(adj.P.Val = p.adjust(p_val, 'fdr')) |>
  filter(gene == 'TRPM2') |>
  mutate(adj.P.Val, logFC = avg_log2FC,
         dataset = 'FLU_CNP0001102', .keep = 'none') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/FLU_CNP0001102.csv')

m2_covvhc <- sobj_bgi |>
  FindMarkers(group.by = 'group', ident.1 = 'COV', ident.2 = 'Ctrl',
              min.pct = 0.001)

m2_covvhc |>
  as_tibble(rownames = 'gene') |>
  mutate(adj.P.Val = p.adjust(p_val, 'fdr')) |>
  filter(gene == 'TRPM2') |>
  mutate(adj.P.Val, logFC = avg_log2FC,
         dataset = 'COVID19_CNP0001102', .keep = 'none') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/COVID19_CNP0001102.csv')

ct_m2_fluvhc <- sobj_bgi |>
  FindMarkersAcrossVar(split.by = 'cell_type', group.by = 'group',
                       ident.1 = 'Flu', ident.2 = 'Ctrl',
                       features = 'TRPM2', min.pct = 0)

ct_m2_fluvhc |> View()

ct_m2_covvhc <- sobj_bgi |>
  FindMarkersAcrossVar(split.by = 'cell_type', group.by = 'group',
                       ident.1 = 'COV', ident.2 = 'Ctrl',
                       features = 'TRPM2', min.pct = 0)

ct_m2_covvhc |> View()

sobj_bgi |>
  DotPlot2d('TRPM2', Stage, cell_type)

## Pseudobulk ---------
pbulk_bgi <- sobj_bgi |>
  AggregateExpression(return.seurat = T,
                      group.by = c('orig.ident','group'))

pbulk_bgi_res <- pbulk_bgi |>
  FindMarkers(group.by = 'group', ident.1 = 'COV', ident.2 = 'Ctrl',
              test.use = 'DESeq2')

pbulk_bgi_res |>
  as_tibble(rownames = 'gene') |>
  mutate(fdr = p.adjust(p_val, 'fdr')) |>
  filter(gene == 'TRPM2')

# blish ------------
sobj_blish <- read_rds('DE_cells/data/blish_cleaned.rds') |>
  UpdateSeuratObject()

sobj_blish$Status |> unique()

sobj_blish |>
  DotPlot2d('TRPM2', Status, cell.type.coarse)

m2_blish_covvhc <- sobj_blish |>
  FindMarkers(group.by = 'Status', ident.1 = 'COVID')

m2_blish_covvhc |>
  as_tibble(rownames = 'gene') |>
  mutate(adj.P.Val = p.adjust(p_val, method = 'fdr')) |>
  filter(gene == 'TRPM2') |>
  mutate(adj.P.Val, logFC = avg_log2FC,
         dataset = 'COVID19_GSE174072', .keep = 'none') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/COVID19_GSE174072.csv')

m2_blish_celltype_covvhc <- sobj_blish |>
  FindMarkersAcrossVar(split.by = 'cell.type.coarse', group.by = 'Status',
                       ident.1 = 'COVID', features = 'TRPM2')

m2_blish_covvhc |>
  mutate(adj.P.Val = p_val*26361, logFC = avg_log2FC,
         dataset = 'COVID19_GSE174072', .keep = 'none')

m2_blish_celltype_covvhc |> DT::datatable()

## Pseudobulk ---------
pbulk_blish <- sobj_blish |>
  AggregateExpression(return.seurat = T,
                      group.by = c('orig.ident','Status'))

pbulk_blish_res <- pbulk_blish |>
  FindMarkers(group.by = 'Status', ident.1 = 'COVID', test.use = 'DESeq2')

pbulk_blish_res |>
  as_tibble(rownames = 'gene') |>
  mutate(fdr = p.adjust(p_val, 'bonferroni'))

# GSE149689 (COV+FLU scRNA) -------
getGEOSuppFiles('GSE149689', makeDirectory = F, fetch_files = F)

sobj <-
list.files('mission/SLE_TRPM2_MfMo/data', 'GSE149689', full.names = T) |>
  read_geo_supp('GSE\\d+')

covflu_meta <-
read_delim('mission/SLE_TRPM2_MfMo/data/gse149689_meta.txt', delim = ' - ',
           col_names = c('suffix','id','sample'), quote = "'") |>
  mutate(group = str_extract(sample, 'nCoV|Flu|Normal'), id = NULL,
         sample = str_remove(sample, ' scRNA-seq'))

sobj <- sobj |>
  mutate(suffix = str_remove(.cell, '.+[A-Z]+') |> as.double()) |>
  left_join(covflu_meta) |>
  mutate(orig.ident = sample) |>
  PercentageFeatureSet('^MT-', col.name = 'mito_ratio')

Idents(sobj) <- sobj$orig.ident

sobj |> VlnPlot('mito_ratio', pt.size = 0)

sobj <- sobj |> filter(mito_ratio < 10)

sobj <- sobj |> quick_process_seurat()

flu_m2 <- sobj |>
  filter(seurat_clusters %in% c(17,1,4,11,16,9)) |>
  FindMarkers(group.by = 'group', ident.1 = 'Flu', ident.2 = 'Normal',
              features = 'TRPM2')

flu_m2 |>
  as_tibble(rownames = 'gene') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE149689_flu_m2.csv')

sobj |>
  DotPlot(pbmc_markers, cols = 'RdBu', cluster.idents = T)

ncov_m2 <- sobj |>
  filter(seurat_clusters %in% c(17,1,4,11,16,9)) |>
  FindMarkers(group.by = 'group', ident.1 = 'flu', ident.2 = 'Normal',
              features = 'TRPM2')

ncov_m2

# GSE101702 (FLU array) ---------
schughart19 <- pluck_geo('GSE101702')

schughart19 |> fData() |>
  as_tibble()

schughart19$group <- ifelse(schughart19$severity.ch1 == 'hlty_ctrl',
                            'ctrl', 'flu')

s19_ctrl_svre <- schughart19[,schughart19$severity.ch1 != 'flu_mod']

s19_ctrl_svre$group |> unique()

s19_ctrl_mod <- schughart19[,schughart19$severity.ch1 != 'flu_svre']

s19_deg <- schughart19 |>
  geo_limma(gene_col = 'GENE_SYMBOL')

s19_deg |>
  filter(gene == 'TRPM2')

s19_svre_ctrl_deg <- s19_ctrl_svre |>
  geo_limma(gene_col = 'GENE_SYMBOL')

s19_svre_ctrl_deg |>
  filter(gene == 'TRPM2') |>
  mutate(logFC = -logFC) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE101702_flu_m2.csv')

s19_mod_deg <- s19_ctrl_mod |>
  geo_limma(gene_col = 'GENE_SYMBOL')

s19_mod_deg |>
  filter(gene == 'TRPM2') |>
  mutate(logFC = -logFC) |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE101702_flu_m2.csv')

# GSE111368 (FLU array) ----------
dunning18 <-
  pluck_geo(file = 'mission/SLE_TRPM2_MfMo/data/GSE111368_series_matrix.txt.gz')

dunning18 |> fData() |>
  as_tibble()

dunning18 |> pData() |>
  glimpse()

dunning18$flu_type.ch1 |> table()
dunning18$t1severity.ch1 |> table()

dunning18$group <- ifelse(dunning18$flu_type.ch1 == 'HC', 'HC', 'Flu')

d18_severe1 <- dunning18[, dunning18$t1severity.ch1 %in% c('HC','1')]

d18_deg <- dunning18 |> geo_limma()

d18_deg |>
  filter(gene == 'TRPM2') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE111368_flu_m2.csv')

d18_severe1_deg <- d18_severe1 |> geo_limma(use_vooma = F)

d18_severe1_deg |>
  filter(gene == 'TRPM2') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE111368_flu_m2.csv')

# GSE157103 (COV RNAseq) ----------
jaitovich20 <- download_ncbi_counts('GSE157103')

jaitovich20

jaitovich20_meta <- download_ncbi_meta('GSE157103')

j20_tdb <- jaitovich20 |>
  pivot_longer(-1, names_to = 'geo_accession') |>
  left_join(jaitovich20_meta) |>
  mutate(group = str_extract(title, '.+VID'),
         severity = str_extract(title, '(Non|)ICU'),
         sex = str_extract(title, '(fe|)male')) |>
  tidybulk(.sample = title, .transcript = Symbol, .abundance = value)

j20_tdb <- j20_tdb |> quick_process_bulk(group = group, skip_scale = T)

j20_qlf <- j20_tdb |>
  test_differential_abundance(~ 0 + group, contrast = 'groupNONCOVID-groupCOVID',
                              omit_contrast_in_colnames = T)

j20_m2 <- j20_qlf |>
  pivot_transcript() |>
  filter(Symbol == 'TRPM2') |>
  mutate(adj.P.Val = FDR, dataset = 'COVID19_GSE157103') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/GSE157103_cov_m2.csv')
